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Abstract 



In a recent work [8] a numerical method has been proposed to simulate off- 
0> I equilibrium zero-temperature parallel dynamics for the SK model without 

^P ■ finite size eff'ects. We study the extension of the method to non-zero tem- 

perature and sequential dynamics, and analyze more carefully the involved 
■^ I computational problems. We find evidence, in the glassy phase, for an al- 

^^ ■ gebraic relaxation of the energy density to its equilibrium value, at least 

at large enough temperatures, and for an algebraic relaxation of the mag- 
netization to zero at non-zero temperatures, with an exponent directly 
proportional to the temperature. 



1. Introduction 

In the last years the study of off-equihbrium dynamics for spin glasses has become of great 
interest, experimentally, numerically and theoretically (as a subjective and non exhaustive 
list of references see, e.g., [1, 2, 3]). The analytical studies are based mostly on the dynamic 
functional integral method of de Dominicis [4] (DFM in the following) by which, in the 
infinite range limit where a mean field solution is exact, the dynamics can be expressed as 
a stochastic equation of motion of a single spin with self consistent nucleus of evolution 
and gaussian noise. In the equilibrium case these equations can be handled imposing time 
translation invariance and a suitable way of recovering partial validity of the fluctuation 
dissipation theorem, which is violated in the spin glass phase because of ergodicity breaking. 
For such an analysis we recall the fundamental work of Sompolinsky and Zippelius on the 
SK model [5] and recent works on the sSG model with p-spin interaction [6]. 

In spite of the great effort of the last years, nowadays the analysis of the off-equilibrium 
dynamics of mean field models has not been developed to the same extension as its equi- 
librium counterpart. From the numerical point of view the situation is reversed. An initial 
state for the dynamical evolution is very difficult to be realized at equilibrium, whereas it 
is easy in the off-equilibrium case: the archetypical off-equilibrium state can be the one 
in which all the spins of the system point in the same direction, or else the one in which 
the spins point independently in random directions. The dynamical equations obtained 
by DFM in the mean field case are single-spin equations at thermodynamical limit: this 
provides, in principle, a powerful tool to simulate by numerical methods off-equilibrium 
dynamics without finite-size effects. However, these equations are not suitable for a nu- 
merical simulation, as they represent the evolution of continuous spins with continuous 
time. A few years ago in [7] it has been developed an approach to the (parallel) dynamics 
of mean field models in the case of Ising spins and discrete times: this approach leads 
to single-spin equations completely analogous to the ones obtained by DFM in the con- 
tinuous case. Recently, for the SK model at zero temperature, in [8] an implementation 
of these single-spin equations, written in the form of stochastic finite- difference equations 
with independent forcing Gaussian noises, has been presented. 

The aim of this paper is to study more carefully the computational problem involved 
in the implementation, introducing a thermal noise at finite temperature, and discussing 
the application of the method to sequential dynamics too. We prove that the method is 
reliable for studying sequential dynamics at every temperature, as far as the relaxation 
of physical quantities of the SK model is concerned. As a first result, in the simplest 
experimental situation of an evolution at constant temperature and zero magnetic field, 
starting from an initial quench from the high-temperature paramagnetic state, we study 
the relaxation of the energy density of the model, which was ignored in [7, 8], together 
with the magnetization, in the whole range of temperature covering the glassy phase. We 
find evidence for an algebraic relaxation of the energy density to its equilibrium value, at 
least at large enough temperatures, and for an algebraic relaxation of the magnetization to 
zero at non-zero temperatures, with an exponent directly proportional to the temperature. 



This behaviour for the magnetization is weU estabhshed experimentaUy on real spin glasses 
(see, for instance, [14, 15]), but it has never been previously tested numerically for mean 
field models with sufficient reliability, since so far in all the numerical simulations it has 
been difficult to take finite-size effects into account (see, for instance, [2]). 

The performance of the algorithm, in terms of CPU-time and memory requirements, 
is not so exceptional to qualify it as the panacea of the dynamical simulations for the 
SK model. However we believe that it is not a remote possibility to simulate situations 
directly comparable with the experimental ones [1] , which present more delicate schedules 
of variation of the physical parameters (e.g. temperature and applied magnetic field) and 
make observable a phenomenology not confined to a simple relaxation of the physical 
quantities. 

The outline of the paper is the following. In section 2 we discuss the generalization 
of the dynamical equations of [7, 8] to finite temperature, and we justify their use also for 
a sequential dynamics, in section 3 we enter in the details of implementation and precise 
the involved computational problems, in section 4 we present the results of our runs. 



2. Dynamical equations and physical quantities of interest 

The Sherrington-Kirkpatrick model for spin glasses ([9], SK model in the following) is a 
model for Ising spins, interacting via a disordered mean field Hamiltonian 

HjW] = --^a,J^^jaj (1) 

Here the (symmetric) couplings Jij are Gaussian independent random variables, i.e. 

J^,J = JJ,^ , rf/x(J,,,) = y^e-^^-/2rfJ,,, (2) 

{N is the size of the sample; the variance of the J's is 1/A^ in order to have finite Hamilto- 
nian density). We denote by ( . ) the thermal average with respect to the Hamiltonian (1) 
at fixed J, and by ( . ) the extra average over the couplings with the measure (2). 

Our aim is to study a suitable dynamics for this model. The interesting physical 
quantities will be densities of time-varying extensive quantities. In order to avoid misun- 
derstandings we recall that the couplings J are constant along the dynamical evolution. 
Typically we shall consider the spin-spin correlation function 



1 ^ 

C(t,t')= lim ^$](a.(t)a,(tO) 



i=l 



and the magnetization 



-w = JL-^^E(-^W) 



i=l 



If ai(t = 0) = 1 for each site i it is 

m{t) =C{t,0) 

The energy density is obviously given by 
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e(t)= lim -{Hj[am 

We recall now the results of [7, 8]. At zero temperature a parallel dynamics for the 
model is defined by the equation 

(T,(t + 1) = sign[^ AjCTjit)] Vz (3) 

Averaging over the J's in the limit A^ -^ oo, from equation (3) in [7, 8] the following 
stochastic single-spin equation is derived: 

a{t+l) =sign[h{t)] 

h{t) = J2^t,sVi^)+J2^t^^^('^ ^^^ 

s<t s<t 

where ry is a temporary uncorrelated normal noise, i.e. the ?7(s)'s at every time s are 
independent random variables, each of which is Gaussian distributed with zero mean and 
variance one. The nuclei A and K in (4) will be fixed by suitable self-consistency equations 
(see (5) below). 

The equation (4) is equivalent to the (3) as far as mean values of density of extensive 
quantities are concerned, in a meaning that will be clarified in the following. In principle 
equation (4) can be solved for each realization of the noise ?7, thus giving the time evolution 
a{t) in function of the noise history r]{s). We denote by ( . )^ the average over the noise 
T]. The nuclei A and K are fixed by the self-consistency equations 

At,t' = for t' > t , Yl ^t,sAt',s = {a{t)a{t'))^ (5.1) 

s<min(t,t') 

and 

Kt,t' = hr t' > t , Yl Kt,sAs,t' = Ht)r]{t'))r, (5.2) 

t'<s<t 

Translated in the formalism of the single-spin equation (4) the spin-spin correlation func- 
tion and the magnetization become 

C{t,t') = {a{t)a{t')), 

m{t) = {a{t)), ^ ^ 

Again, if a (t = 0) = 1, it is 

m(t) =C(t,0) 
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A very high precision plot of the relaxation of m{t) versus t until t = 100 with initial 
condition m(0) = 1, is shown in [8]. The data in [8] can well be fitted with an algebraic 
relaxation of the type 

m{t) ~ nioo + At'"" 

and the quoted results for moo and a are 

moo = 0.184 ±0.002 a = 0.474 ± 0.005 
The result for nioo agrees quite well with what is declared by [7] 

moo = 0.23 ±0.02 
It is commonly believed that at T 7^ this relaxation form becomes 

m(t) ~ cr^ 

with zero remanent magnetization and 6 ccT. We shall come back to this point in section 4. 
As far as the energy density is concerned, a direct translation in the formalism of the 
single-spin equation (4) would result in 

e{t) = -l{a{t)h{t)), (7) 

but this expression is identically zero. In fact, a peculiar characteristic of the dynamics (4), 
already pointed out in [7, 8], is that 

{a{t)a{t'))rj = if |t-t'| isodd ,. 

(a(t)?7(t'))r, = if |t-t'| iseven ^^ 

and it is easy to see that it implies that the form (7) for the energy density gives e(t) = 0. 
In [7, 8] this inconsistency problem does not arise, because the energy density is ignored 
at all. We shall see soon that the correct expression is 

e{t) = -lHt)h{t - 1)), = -^{a{t - l)h{t))r, (9) 

The source of the problems is the fact that the parallel dynamics (3) is not a correct 
thermal dynamics, even at T = 0. In other words it does not drive the system to its 
thermal equilibrium, as does its sequential counterpart. The dynamics (3) is of practical 
interest for models of attr actor neural networks, but at a first sight is not such for the SK 
model, considered as a thermodynamical model. Keeping this in mind it is not difficult to 
accept the fact that the energy density is identically constant in time, i.e. the dynamics 
is not dissipative. We can get rid of this unpleasant situation by considering a slightly 
different model, the Little model [10], for which a parallel dynamics too is a good thermal 
dynamics. We shall see soon that the dynamics deffned by (3), which is a parallel (i.e. 
a bad) dynamics for the SK model, is instead a good thermal dynamics for the Little 



model. In order to define the model we double the phase space considering, instead of A^ 
Ising spins {o'i}i, 2N spins in two sets {a^ji and {ri}i, for i = 1..N, interacting via the 
disordered mean field Hamiltonian 

Hj[o-, t] = - ^ (^iJijTj (10) 

Again, the (symmetric) couplings J are Gaussian random variables with the same mea- 
sure (2). Now, it is easy to see that the (semi)parallel dynamics 

a,{t + 1) = sign[^ Jz,jrj{t)] , n{t + 1) = n{t) Vz 

'^^ (11) 

r,(t + 2) = sign[^(T,(t + 1) Jj- ,] , (T,(t + 2) = a,{t + 1) Vi 

is equivalent to a sequential one, thus providing a good thermal dynamics. 

In order to recover equation (3) we change slightly the notation. Starting from t = 0, 
(11) gives Ti{t = 1) = Ti{t = 0), whereas the value of ai{t = 1) is not trivial. Then 
ai{t = 2) = ai{t = 1) and Ti{t = 2) is not trivial, and so on. In general the value of 
the r's at odd times and the value of the a's at even times is trivial, so these variables 
are not interesting, as far as the dynamical description of the system is concerned. Now, 
we rename tout court the remaining r's (the ones at even times) simply by a. Thus we 
are left with an unique set of dynamical variables o"i(t), and, rewritten in terms of it, the 
equation (11) is reduced to the previous (3). 

At this point we are left with the following situation: the (formally) parallel dynam- 
ics (3) is a good thermal dynamics, at zero temperature, for the Little model. We shall 
see soon that the introduction of a finite temperature thermal noise does not present any 
difficulty. However, we must firstly justify why we feel us authorized to study the Little 
model instead of the SK model. In fact, the Little model in itself is of little physical inter- 
est, but it has ben shown [11], both by a replica approach and numerical simulations, that 
in the limit N ^ oo the static properties of the two models are the same, i.e. their ther- 
modynamical behaviours coincide. In particular, as far as the energy density is concerned, 

it is 

e= lim l(iff^)[a]) = l lim l(iff ^«'^)[ct,t]) (12) 

Now, with our notation, the dependence on time of the Little hamiltonian (10) is 
i:fi^^"'^)[a(t),r(t)] = -5]a,(t)J,,,r,(t) = -5^a,(t)J,,,a,(t-l) 

Thus it is easy to see that, in the formalism of the single-spin equation (4), the energy 
density becomes 

eL^ttleit) = -{a{t)h{t - 1))^ = -{a{t - l)h{t))n 

From this expression, and relation (12) between the energies in the Little and SK models, 
we recover the previously anticipated equation (9) for the correct form of the energy density 
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( "correct" in the sense that we can compare correctly the dynamical results with the static 
ones for the SK model). 

To extend the dynamical equations to finite temperature, it is sufficient to introduce in 
the equation (3) a temperature dependent local effective ffeld which simulates the thermal 
agitation. This ffeld, because of the fact that it does not depend on the J's, is carried out 
without modiffcations during the derivation of the single-spin equation (4) starting from 
the (3). The resulting equation, which generalizes the (4) of [8], is 

a{t + 1) = sign[h{t)] 
h{t) = Y, AtMs) + E ^M^(«) - T ■ atMt - 1) (13) 

s<t s<t 

where the nuclei A and K will again be ffxed by the (5), and tj is again a temporary 
uncorrelated normal noise. The thermal induced noise ^ is a new temporary uncorrelated 
noise with law 

P(0 = 0(0 ■ 2e-2^ (14.1) 

which represents a thermal dynamics of Metropolis type, or 

P(0 = ^[l-tanh2(0] (14.2) 

which represents a thermal dynamics of heath-bath type. Obviously the presence of this 
second noise forces us to slightly change the notation in the above relations (5), (6) and (9), 
now denoting with ( . )^ ,e the average over both noises. 

Let us stress the fact that in (13) the thermal ffeld is proportional to a{t — 1) and 
not to o"(t). This is a consequence of the fact that actually we are studying the dynamics 
of the Little model and not of the SK model: when we are updating a spin, deciding if 
we ffip it or not, its previous value is actually cr(t — 1), whereas a{t) is only a short- hand 
notation for T{t). Obviously in the case of the heat-bath dynamics the thermal ffeld can 
be reduced simply to T ■ ^(t), as in this case there is no need of a previous spin to ffip: in 
fact, the distribution law of the heat-bath noise (14.2) is even in ^. For generality's sake we 
maintain the form (13) for the dynamical equations, whereas for practical use it is better 
to pursue in each case simplicity in place of generality. In next section we shall see more 
carefully how to implement a simulation for the single-spin equation (13), and in section 4 
we give the results of our runs. 



3. Details of implementation and computational problems 

To implement a simulation of the stochastic ffnite-difference equation (13) we follow the 
main lines of [8]. We generate randomly Nt different noise trajectories, and we deffne the 
average over the noise simply by 



Nt 

l/rp 



(a(t)a(t')),,c = ]^E^^^^W^^^^(^') 



1 



and by analogous expressions for the other interesting quantities. Here a^^\t) denotes the 
evolution of a in the z-th noise history. Let us stress that the run-time computation of 
certain averages is needed for the evolution itself, since the nuclei A and K are defined in 
terms of averages over the noise. 

We shall be more explicit in a while, let us before anticipate two comments. Firstly, 
in the runs which we performed we kept the temperature fixed, and zero magnetic field. 
Changing T and h during the evolution would have been pointless, as our care was mainly 
devoted to test the reliability of the method. 

Secondly, the choice of the initial conditions for the evolution is not straightforward. 
If in the equation (13) there was not the term a{t — 1), as it is for T = or for finite 
temperature heat-bath dynamics, then we would need only the initial condition at t = 0, 
i.e. {a^^\t = 0)}i=i..Ary. In this case, as pointed out in [7], the choice of the initial condition 
is to a great extent arbitrary. In fact, because of the symmetry of the equation (13) under 
simultaneous inversion of a and both noises ry and ^, it is not difficult to see that the 
average of a product of an even number of fields {a and/or the noises) does not depend on 
the initial conditions. The average of a product of an odd number of fields, on the other 
hand, depends on the initial conditions only through a factor which is simply given by the 
initial magnetization 



Nt 



m(0) = -^ V(7(')(t = 0) 



=1 

From the physical point of view the initial condition for the evolution represents a sudden 
quench from a high temperature paramagnetic state which is in equilibrium in a field that 
sets the value of the magnetization. 

In the general case the presence in (13) of the term a{t — 1) forces us to impose two 
distinct initial conditions: besides {a^'^\t = 0)}i we need also {a^'^\t = — l)}i. At this 
point it is not clear at all neither the exact physical meaning of this additional initial 
condition, nor how to keep a priori under control the dependence of the evolution on the 
initial conditions. However, it is reasonable to argue that the additional initial condition 
at t = — 1 has no real significance. We have therefore fixed the initial conditions in the 
way we believe to be the most natural one: 

a«(t = 0) = l yi^L.Nr (15.1) 

and 

aW(t = -l) = l Vz = 1..A^t/2 

(15 2) 
aW(t = -l) = -l \/i = NT/2 + l..NT 

We have not yet performed any systematic test on the dependence of the evolution on 
the initial conditions. Let us stress, by the way, that with the choice (15.2) the dynamics 
preserves the same property of selecting the parity (8) as in the zero temperature case. 

We shall now sketch with some details the algorithm used in the simulation. Besides 
the already introduced notation, we denote respectively with S and X the spin-spin and 
the spin-noise correlation functions, i.e. 

St,t' = (a(t)(7(t'))^,^ Xt^t' = Ht)r]{t'))r,,^ 



The main steps are the following: 

1.1) fix the values of: 

- type of thermal dynamics (Metropolis or heat-bath) 

- temperature T 

- number of noise histories Nt 

- terminal time point of the evolution tmax 

- initial seed of the random generator 

1.2) for each history i = 1..Nt, fix the initial conditions for a"(0) and a"(— 1) as in (15) 

1.3) for each history i = 1..Nt, generate the t = values of the noises ?7(0) and ^(0) in 
according, respectively, to a normal law and to the law (14) 

1.4) set the t — values of the nuclei and of the correlation functions, i.e. 

'S'oo = 1 ^00 = 1 -'^oo = Kqo = 

1.5) set t = 

2) let t be the current value of the time step; if t > tmax then stop the main execution 
and go to step 8) (output of results) 

3) for each history i = 1..Nt, compute a{t + 1) in according to the equation (13), by 
using the nuclei A and K already computed in the previous steps 

4) compute 

. Nt 

5i+i,, = — Va«(t + l)a«(s) s<t + l 



Ss,t+1 — St + l,s 7 St+l^t + l — 1 



and 



Nt 



1=1 

( Xs,t+i = , Xt_|_i^t_|_i = 



(the correlation functions Ss,s' and Xs^s' for each s,s' < t are already known from the 
previous steps) 
5.1) compute the nucleus A^+i ^ by inverting its defining equation (5.1); this gives the 
equation 

^t+i,t' = 2_^ St+i^sA^,,. t' <t 

s<t' 



At+i,t+i= k-Y.{At+i,t'y 



t'<t 



which is well-posed since we already know from the previous steps Ag^s', for each 
s,s' <t 
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5.2) compute the nucleus i^t+i,s by inverting its defining equation (5.2); this gives the 
equation 

Kt+i,t' = 5Z ^t+i,sA~l, t' < t 

t'<s<t+l 

Kt+i,t+i = 

which again is well-posed since we already know from the previous steps Ag^s', for 
each s,s' <t 

6) for each history i = 1..Nt, generate the values of the noises r]{t + 1) and ^{t + 1), in 
according, respectively, to a normal law and to the law (14) 

7) set t = t + 1 and go to step 2) 

8) output the results of the run: 

- St,t' and Xt,t' for t, t' < tmax 

- m{t) = St,Q 

- e{t) = -1/2 J2s<t At-i,sXt,s + At,sXt-i,s 

Two technical remarks. Firstly, the random generator for the noises is the usual machine 
generator RAN of the VAX-VMS systems. It is useful to generate two distinct random 
series for the two noises r] and ^. The generator RAN gives an uniform distribution between 
and 1, but there is no difficulty to obtain, starting from it, the required distributions 
for ry and ^. Secondly, in steps 5.1 and 5.2 the computation of A~^ is easily made by the 
Gauss method thanks to the fact that A is a (lower) triangular matrix. 

As far as the C PU-time requiring is concerned, we note that the most expensive steps 
are the 3) and the 4), which cost O(A^t^) multiplications and sums, and the 5)'s, which 
cost O(t^). Summing over t until t = tmax we have 



CPUr^c-NTtl^^ + c'-t 



max ' "^ ^max 



Our runs were performed on a DEC-VAX 6000, running VAX Fortran, where for Nt = 1000 
and tmax = 200 the required CPU-time is about 1 minute. Concerning the memory 
occupation, the matrices S, A and X require 0{tmax) memory, whereas K can be stored 
in a scratch array of 0{tmax), since at time step t we need only {Kt,s}s<t- The main 
memory occupation stems however from the complete Nt histories of a and tj, whereas 
^ is needed only at current time step t. This gives an occupation of 0{NTtmax)- In 
conclusion, reducing the memory occupation as much as possible, we have 



Mem/bytes ~ 5 ■ Nxtmax + 3 ■ t 



2 

max 



which gives an occupation of about 1.1 Mbytes for Nt = 1000 and tmax = 200. 

Roughly speaking, the CPU-time sets a limitation for tmax, whereas the memory does 
it for Nt- If we want to reach a higher value of tmax we only need to wait more time. At a 
first glance, the significance of Nt is not so clear. It obviously affects the precision of the 
computation, but the error due to finite Nt is not a purely statistical one, of the order of 
I/^/Nt, because average quantities enter in the evolution itself throughout the nuclei A 
and K. We cannot thus assume a priori that, repeating more times the same run with a 
low value of Nt, we shall obtain a higher precision. 
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It is therefore necessary to understand, at least empirically, the dependence of the 
computation results on A^^. This for two reasons: the inherent limitations on a particular 
machine force us to keep the value of Nt relatively small (Nt = 32, 000 at maximum, for 
our runs); moreover, even if we would dispose of a bigger machine, we ought to control a 
priori the reliability of a computation, in order to be confident of its results. Thanks to 
the availability of a Cray, in fact, the authors of [8] could keep Nt = 1,000,000. Such 
a value appears sufficiently high, although for no a priori reason): we shall come back to 
this point in the next section. By the way, though it is not explicitly declared, probably 
in [8] a slightly different implementation was used. In fact, Nt — 1, 000, 000 would impose 
a memory occupation of about 1 Gbyte, which is enormous even for a Cray. To reduce 
the memory occupation one could simply avoid to store in memory the complete a and 
?7 histories: one could maintain an array of 0{tmax) with the initial conditions for a and 
the initial seeds of the random generator for ry, and at every time step t one computes, or 
generates, ex novo the a and tj histories until that time, by using the nuclei A and K already 
computed and stored in memory. This makes the memory occupation independent on the 
value of Nt, thus releasing any limitations on it, but increases the CPU-time requirement 
by a factor 0{tmax)- On a VAX this would bring to about 25 days the amount of CPU- 
time spent in a run with Nt = 1,000,000 and tmax = 100, but a Cray Y-MP could be 
even 100 times faster than a VAX, thus reducing the time to about 6 hours, which is what 
declared in [8]. 

We shall now anticipate what we regard as the main source of damage linked to 

low values of Nt- Let us consider a given history a^^' . At each time step t there is a 

probability of spin fiip which, in principle, can depend on the whole history {a^^'{s)}s<t- 

Let us suppose, however, that this probability is even under global spin inversion. In this 

case it is easy to see that 

dm / N / N 

_ = -m{t)p{t) 

where p{t) is the marginal spin flip probability at time t. Now, the probability that none 
of the Nt histories will flip at time t is appreciable as soon as p{t)NT < 1. At T = it 
is easy to see that, if at a certain time i it is a^'^'{i -\- 1) = a^'^\i) for each history i, then 
it is a"*-*)(t) = a^^'{i) for each t > i and each history i, and the system freezes. At T 7^ 0, 
because of thermal fluctuations, this freezing can never happen, but we believe that a kind 
of partial freezing forces the system to behave anomalously. In next section we shall show 
such an anomalous behaviour, though we can not prove that this is due to a mechanism 
of the type already described. For the time being, let us note that p{t) is a function which 
decreases to zero as t goes to inflnity, at least if the relaxation of m{t) is slower than 
exponential (an algebraic law or a stretched exponential one), as it is commonly believed. 
Therefore it becomes p{t)NT < 1 as soon as 

t > tcrzt{NT) 

and tcrit is an increasing function of Nt- It is appealing to suppose that 

m{t) ~ t"" 
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i.e. the relaxation of m{t) is algebraic to zero, so that p{t) ~ 1/t and 

tcr^t{NT) ~ Nt (16.1) 

Probably there are some additional complications. In fact at T = 0, as we anticipated, it 

is [7, 8] 

m{t) ~ moo + At~" 

i.e. the algebraic relaxation of m{t) goes to a non-zero value, so that p{t) ~ l/t^"*"" and 

tcr^t{NT) ~ AT^ < A^T (16.2) 

At T 7^ it is expected a relaxation to zero, but even if a kind of freezing mechanism 
works, it is not clear how to implement it. We shall see in the next section that, at least 
qualitatively, the numerical data give some evidence of the presence of a kind of "freezing 
time" which grows with Nt, though it is difficult to distinguish between a situation like 
the one of equation (16.1), where the freezing time grows quite fast with Nt, and the one 
of (16.2). 

4. Numerical results 

We made two series of runs. The first one was devoted to analyze the dependence of the 
results of a run on the value of Nt, in the case T = 0.4 and thermal dynamics of Metropolis 
type, while in the second one we measured the relaxation of the magnetization and of the 
energy density in the whole range of temperatures between 0.1 and 1.5, for both Metropolis 
and heath-bath dynamics. In both series the maximum time reached was tmax = 200, and 
the external magnetic field was kept equal to zero. We tried also to measure the relaxation 
at T = 0, but for the values of Nt we used the system freezes before reaching tmax- 

We discuss firstly the results of the first series. The choice of the Metropolis dynamics 
has no particular motivation. The value of the temperature has been chosen to be not 
too high, in order to have relevant effects of the breaking of the replica symmetry (we 
remember that with our choice of the normalization of the coupling (2) the glassy critical 
temperature is Tg = 1), and not too low, in order to avoid the very slow relaxation and 
the freezing effects of low temperatures. We made several runs for A^^ ranging from 1000 
to 32, 000. Apart from the systematic error due to Nt, which we would investigate, there 
is a natural statistical error of order l/y/Nr, thus for each value of Nt we made many 
independent runs, in order to have comparable statistical errors. More precisely 

Nt = 1000 num. of indep. runs = 80 

2000 40 

4000 20 

8000 10 

16,000 10 

32, 000 10 
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Figure 1. Plot of the relaxation of the energy density — 2e(i) for Metropolis dynamics at temperature 
T = 0.4 and zero magnetic field. The horizontal dotted line represents the equilibrium value. Each curve 
corresponds to a different value of Nt, as it is shown in the inset, and it has been obtained by averaging 
the results of many independent runs, in order to have error bars of comparable magnitude. 

In figure 1 we plot the energy density — 2e(t) versus t for all the values of Nt- The time 
scale from to 200 was divided in intervals of size proportionally scaled by a fixed factor 
(we have chosen 1.5 to have about ten intervals in the considered range), then for each run 
we averaged the data into each interval. In fact, the original data are strongly dependent 
at nearby times, so the form of the plot at small time scales has no significance. With 
our choice the reduced points are equally spaced in logarithmic time scale. This choice 
for the procedure of decimation of the data is due to the fact that we expect an algebraic 
relaxation, i.e. a linear dependence of ln(e — Coo) on Int, and such a functional form is 
invariant under such a transformation, up to an innocent logarithmic translation in the 
time scale. Finally, we made the arithmetic mean of the different runs, and we kept as 
error bars for each datum their standard deviation (Icr). 

The horizontal dotted line is plotted as a reference and represents the equilibrium 
value for the energy density, computed by numerically solving the replica equations for the 
SK model [12]. From figure 1 we see that at a certain time the plot of — 2e(t) starts to grow 
almost linearly and rapidly exceeds its equilibrium value. This anomalous effect is greatly 
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reduced by increasing A^^. It is useless to show a plot of the magnetization analogous to 
the one in figure 1, since in this case the various relaxation curves at different values of 
Nt do not differ sufficiently to make the effect of Nt directly observable. Anyway, if we 
consider the data at fixed t versus Nt, a linear dependence on 1/Nt clearly appears. In 
fact, the data are very well fitted by 

y^,(t)^y^(t) + ^ (17) 

where Y represents indifferently the energy density — 2e or the magnetization m. The 
curves Yoo{t), i.e. the extrapolations to Nt = oo, will be the subject of the next analysis. 
Thanks to the very simple dependence (17) on Nt, to obtain a fairly reliable extrapolation 
of the data Nt = oo, it is sufficient to produce two series of measures at two different 
values of Nt. In the following we shall choose Nt = 8000 and 16,000. The factor A(t) 
in (17) depends roughly linearly on t, more precisely: 

Aenergy{t) ~ 0.6 t 
Amagn.{t) ~ 0.2 t 

In conclusion, the systematic error due to finite Nt, being 0{t/NT), is fairly well manage- 
able: at Nt fixed we can be confident in the results of a run until t remains smaller than 
a (not too small) fraction of Nt. It would be very interesting if we understood why the 
relaxation of the physical observables presents, for finite Nt, such an anomalous behaviour 
as in (17), particularly evident in figure 1 for the energy density. We believe that the un- 
derlying reason is a kind of freezing mechanism as it is sketched at the end of previous 
section, but we could not identify it. 

In figure 2 we show the contour lines of the correlation matrix St,t' in the plane t-t', 
for the different values of Nt, measured in a randomly chosen run of each set. Obviously at 
low values of Nt the plots are very noisy, because of the great statistical errors. Apart from 
this noise, the most evident feature of the plots, which gradually vanishes by increasing 
Nt, is the following. Let us consider the half plane t > t' , being the figure symmetric. 
As soon as t is greater than a certain (not well identifiable) critical time, for each t' , St,t' 
versus t appears to freeze to a value which is almost constant with t. This behaviour is 
strongly reminiscent of the above mentioned freezing. 

The main series of measures was devoted to explore systematically the relaxation 
curves of the energy density and the magnetization in the whole range of temperatures 
between T = 0.1 and T = 1.5. We recall that the glassy transition happens at Tg = 1. For 
each temperature value, and both dynamics (Metropolis or heat-bath) we performed two 
sets of independent runs at two different values of Nt, to extrapolate the data at Nt = oo. 
More precisely 

Nt = 8000 num. of indep. runs = 20 

16,000 10 

(with the exception of T = 0.4 with Metropolis dynamics, where we have already enough 
data). Now, we have pre-analyzed the data as follows. Firstly, for each run, we performed 
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Nt = 4000 



Nt = 8000 
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16, 000 



A^7 



32, 000 
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on the data the previously described procedure of decimation, to obtain a reduced set of 
almost-independent time points. Then, for each value of Nt, we averaged the data from 
the different runs, and the error bars were computed by the mean standard deviation (la). 
Secondly, we extrapolated to Nt — C)0 assuming a linear dependence on 1/Nt as in (17), 
that is we computed, for each time point, 

^cx) = 2 ■ yi6,000 — ^8000 

where the error bars are the obvious ones. For uniformity's sake we performed such an 
extrapolation even for T = 0.4 with Metropolis dynamics, instead of using the result of 
the previously described fitting procedure based on 6 different values of Nt- 

At this point we are left, for each temperature value and dynamics type, with a unique 
set of data giving the time relaxation of the magnetization and the energy density of the 
system. The data are very well fitted by the following functional forms, for T < Tg = 1: 
as expected, the magnetization shows an algebraic relaxation to zero 

m{t) ~ Ct-^ (18.1) 

and the energy density shows an algebraic relaxation to a well-defined non-zero value 

e{t) ^ Coo + C't-^' (18.2) 

A few technical remarks. The fits were performed minimizing the (naturally defined) chi- 
square function. The one for the magnetization can be easily linearized and performed 
by standard methods, the one for the energy density, being highly non-linear, requires a 
package for non-quadratic function minimization and error analysis, in our case the package 
MINUIT of the CERN libraries. To test the reliability of the fits we impose a cutoff tmin 
at low times, we consider the relaxation only for t > tmin, and we check if there is a 
systematic dependence of the fitted parameters on tmin- The fits of the energy density show 
no evident dependence. The ones for the magnetization, for low temperature (T < 0.5), 
show a systematic decrease of the fitted exponent d- This fact is commonly considered a 
sign of the presence of non-negligible corrections to the asymptotic behaviour (18.1). Thus 
we adopt, to take into account such (unknown) corrections, the Berretti-Sokal procedure, 
described in [13, sections 4.2 and 5.3]: we fit the data with the functional form 

m{t) =Ct-^{l + -) 

V 

for various Gxed values of B, and we look for the range in B for which there is no systematic 
dependence of the fitted parameters d and C on tmin- In this range of S, for each fitted 
parameter we select the maximum and the minimum values: the best fit for the parameter 
will be simply the arithmetic mean of these values, whereas their difference will be con- 
sidered as the systematic error due to unconsidered corrections to the leading behaviour, 
or to imperfect knowledge of the form of the corrections (95% subjective confidence limit 
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Figure 3. Plot of the asymptotic value — 2eoo of the energy density versus T for Metropolis dynamics at 
zero magnetic field. The data where obtained by a fit with the algebraic relaxation form (18.2) for T < Tg, 
whereas in the paramagnetic phase for T > Tg the energy density relaxes almost instantaneously to its 
asymptotic value. The full line is the plot of the energy density for the SK model at equilibrium [12]; the 
tail for T > Tg has the usual paramagnetic form e = — 1/2T. The error bars from the fit are negligible 
with respect to the dimension of the dots. 



as defined in [13, footnote 17]), in addition to the usual statistical error for the fit, at yo/o 
confidence level (2a). 

Now, the results of the whole analysis can be summarized, for both dynamics, by 
three functions of temperature: the exponents 5 and 5' of the algebraic relaxation of the 
magnetization and the energy density, and the asymptotic value Coo of the energy density. 
The interesting results are the following. 

- The energy density relaxes to its equilibrium value, at least for high enough tempera- 
tures (T > 0.2 -^ 0.3). The plot of — 2eoo versus T, for Metropolis dynamics, is shown 
in figure 3; the corresponding one for heat-bath dynamics is very similar. The full 
line is the plot of the energy density for the SK model at equilibrium, obtained by 
an exact numerical solution of the replica equations [12]. Recall that for T > T^ an 
exponential relaxation is expected, instead of the algebraic one of equation (18.2), 
and this can not be observed, since the procedure of decimation we followed to reduce 
the data does not leave invariant an exponential relaxation. However, the asymptotic 
value is almost instantaneously reached, so there is no need for a fit. Let us remark 
that, in fact, our estimate for the asymptotic value of the relaxation is slightly lower 
than the equilibrium value (in figure 3 minus two times the energy is plotted, so the 
experimental points appear above the equilibrium curve). This is very probably due 
to some inaccuracy in the extrapolation to Nt = oo. 

- The exponent 5' of the relaxation of the energy density is almost constant with T. 
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Figure 4. Plot of the exponent 5' of the algebraic relaxation of the energy density versus T, for both 
Metropolis (■) and heat-bath (D) dynamics at zero magnetic field. The data where obtained by a fit with 
the form (18.2) for T < Tg, whereas in the paramagnetic region the relaxation is exponential. The two 
series of data hereby displayed are slightly shifted for mere convenience of representation. The error bars 
from the fit are at 95% confidence level (2a). 

A plot of 6' versus T, for both dynamics, is shown in figure 4; let us remark that 
for Metropolis dynamics it is 5' ~ 1, slowly increasing with T, whereas for heat-bath 
dynamics it is 6' ~ 0.7, again slowly increasing with T, but with a jump at T = T^ = 1 
where it is 5' ~ 1. This behaviour of the exponent 6' is somewhat anomalous, since 
it differs from one dynamics to the other, whereas we expect an universal behaviour. 
Moreover, though our results agree qualitatively with previous simulations [2] , it would 
be nice if the exponent of the relaxation go to zero with T. We shall come back to 
this point in a while. 
- The exponent 6 of the relaxation of the magnetization is proportional to T. This can 
be seen from a plot of 5 versus T, which is shown, for both dynamics, in figure 5. 
Moreover, the data for 6 are very well compatible with a dependence on T of the 
simple form 

5{T) = T/Tg 

though we have not scrupulously tested this by a fit. Noticeably, this feature has been 
observed in experiments on real spin glasses since a long time (see, for instance, [14, 
15]). Strong indications for the validity of this fact also for the mean-field SK model 
come from recent numerical simulations [2], but it has never been previously tested 
numerically with sufficient reliability, since so far in all the numerical simulations it 
has been difficult to take finite-size effects into account. Let us stress that the exact 
numerical solution of the SK dynamics at equilibrium gives an exponent u for the 
algebraic relaxation of the magnetization which tends to a non-zero value when T 
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Figure 5. Plot of the exponent 6 of the algebraic relaxation of the magnetization versus T, for both 
Metropolis (a) and heat-bath (6) dynamics at zero magnetic field. The data where obtained by a fit with 
the form (18.1) for T <Tg, whereas in the paramagnetic region the relaxation is exponential. The dashed 
line 5(T) = T is drawn as a guide for the eye. The error bars are the systematic ones of the Berretti-Sokal 
procedure for T < 0.5 (95% subjective confidence limit as defined in [13, footnote 17]), whereas for T > 0.5 
they are the usual statistical ones from the fit, at 95% confidence level (2a). 

goes to zero [12]. This implies that what is observed cannot be explained by the 
equilibrium features of the SK dynamics. 
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As a concluding remark, a few words of caution are necessary. Though from our data we do 
not observe a strong evidence for the presence of subleading terms in the relaxation form 
of the energy density (18.2), nevertheless there are some anomalies in the behaviour of the 
energy that could be explained by the underlying presence of such subleading corrections. 
In fact, we have already pointed out that at low temperatures (T < 0.3) our estimate of the 
asymptotic value of the energy density is definitely higher than the true equilibrium value. 
Our data (see figure 3) show a clear trend to the value — 2eoo — 1.4 at zero temperature, 
in agreement with what was observed in previous simulations [2]. It is expected that at 
T y^ this value relaxes again very slowly to the true equilibrium value, with an exponent 
which goes to zero with T. However, due to the fact that this further relaxation at finite 
temperature is only a few percent of the main relaxation at zero temperature, one could 
argue that at low temperatures and not sufficiently long times the observed exponent of 
the relaxation is an effective exponent, with three consequences: the effective exponent 
is not universal, and it is different in Metropolis and heat-bath dynamics; it is almost 
constant with T, and remains different from zero even at low temperatures (in [2] it was 
observed 6' ~ 0.6 -^ 0.8 for T ranging from 0.2 to 0.4); we miss the estimate for the 
asymptotic value, which at low temperatures remains higher than the true equilibrium 
value. All those troubles are not present in the relaxation of the magnetization, since its 
zero-temperature asymptotic value is moo — 0.2 [7, 8], so its further finite temperature 
relaxation is much more impressive and much more easy to observe. 
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